Exercise 3: electron vs gamma separation¶
1. Introduction¶
Fig. 13 Difference between an electron vs photon is at the start of the electromagnetic shower, where the photon has a gap. From Wouter Van De Pontseele, ICHEP 2020.¶
Electrons are visible in a LArTPC detector because of the electromagnetic showers that they trigger.
Photons, on the other hand, are neutral (no charge) and thus remain invisible to the LArTPC eyes until they convert into electrons (pair production) or Compton scatter. In both cases, the visible outcome will be an electromagnetic shower.
How can we differentiate the two, then? The answer is in the very beginning of the EM shower. For an electron, this shower will be topologically connected to the interaction vertex where the electron was produced. For a photon, there will be a gap (equal to the photon travel path) until the EM shower start (when the photon becomes indirectly visible through pair production or Compton scatter). That seems simple enough, right? Wrong, of course.
Energetic photons could interact at a distance short enough from the interaction vertex, that we would not be able to see the gap. Or, the hadronic activity might be invisible, because it includes neutral particles or because the particles are too low energy to be seen. In that case the interaction vertex might be hard to identify, and the notion of a gap goes away too. For such cases, fortunately, there is another way to tell electrons from gamma showers. Another major difference is in the energy loss rate at the start of the EM shower. An electron would leave ionization corresponding to a single ionizing particle, whereas a pair of electron + positron coming from a photon pair production would add up to two ionizing particle. Thus, we expect the dE/dx at the beginning of the shower to be roughly twice larger in the case of a gamma-induced shower compared to an electron-induced shower.
Fig. 14 Example from MicroBooNE. Left is the shower \(dE/dx\), right is the gap between the vertex and shower start. From Wouter Van De Pontseele, ICHEP 2020.¶
Why do we care? The difference becomes significant if, for example, you are looking for electron neutrinos. One of the key signatures you would be looking for are electrons.
In this exercise, we will focus on finding the start of EM showers and computing the reconstructed dQ/dx in these segments. Optionally, you could compare that to the result of using automatic PID as predicted by the chain.
2. Setup¶
a. Software and data directory¶
import os, sys
SOFTWARE_DIR = '%s/lartpc_mlreco3d' % os.environ.get('HOME')
DATA_DIR = os.environ.get('DATA_DIR')
# Set software directory
sys.path.append(SOFTWARE_DIR)
b. Numpy, Matplotlib, and Plotly for Visualization and data handling.¶
import numpy as np
import matplotlib.pyplot as plt
import seaborn
seaborn.set(rc={
'figure.figsize':(15, 10),
})
seaborn.set_context('talk')
import plotly
import plotly.graph_objs as go
from plotly.subplots import make_subplots
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=False)
c. MLRECO specific imports for model loading and configuration setup¶
from mlreco.main_funcs import process_config, prepare
import warnings, yaml
warnings.filterwarnings('ignore')
cfg = yaml.load(open('%s/inference.cfg' % DATA_DIR, 'r').read().replace('DATA_DIR', DATA_DIR),Loader=yaml.Loader)
process_config(cfg, verbose=False)
/usr/local/lib/python3.8/dist-packages/MinkowskiEngine/__init__.py:36: UserWarning:
The environment variable `OMP_NUM_THREADS` not set. MinkowskiEngine will automatically set `OMP_NUM_THREADS=16`. If you want to set `OMP_NUM_THREADS` manually, please export it on the command line before running a python script. e.g. `export OMP_NUM_THREADS=12; python your_program.py`. It is recommended to set it below 24.
Config processed at: Linux tur015 3.10.0-1160.42.2.el7.x86_64 #1 SMP Tue Sep 7 14:49:57 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux
$CUDA_VISIBLE_DEVICES="0"
d. Initialize and load weights to model using Trainer.¶
# prepare function configures necessary "handlers"
hs = prepare(cfg)
dataset = hs.data_io_iter
Welcome to JupyROOT 6.22/09
Loading file: /sdf/home/l/ldomine/lartpc_mlreco3d_tutorials/book/data/mpvmpr_012022_test_small.root
Loading tree sparse3d_reco
Loading tree sparse3d_reco_chi2
Loading tree sparse3d_pcluster_semantics_ghost
Loading tree cluster3d_pcluster
Loading tree particle_pcluster
Loading tree particle_mpv
Loading tree sparse3d_pcluster_semantics
Loading tree sparse3d_pcluster
Loading tree particle_corrected
Warning in <TClass::Init>: no dictionary for class larcv::EventNeutrino is available
Warning in <TClass::Init>: no dictionary for class larcv::NeutrinoSet is available
Warning in <TClass::Init>: no dictionary for class larcv::Neutrino is available
Since one of the GNNs are turned on, process_fragments is turned ON.
Fragment processing is turned ON. When training CNN models from
scratch, we recommend turning fragment processing OFF as without
reliable segmentation and/or cnn clustering outputs this could take
prohibitively large training iterations.
Since one of the GNNs are turned on, process_fragments is turned ON.
Fragment processing is turned ON. When training CNN models from
scratch, we recommend turning fragment processing OFF as without
reliable segmentation and/or cnn clustering outputs this could take
prohibitively large training iterations.
Ghost Masking is enabled for UResNet Segmentation
Ghost Masking is enabled for MinkPPN.
Restoring weights for from /sdf/home/l/ldomine/lartpc_mlreco3d_tutorials/book/data/weights_full_mpvmpr_012022.ckpt...
Done.
Let’s load one iteration worth of data into our notebook:
data, result = hs.trainer.forward(dataset)
Segmentation Accuracy: 0.9824
PPN Accuracy: 0.8162
Clustering Accuracy: 0.9279
Shower fragment clustering accuracy: 0.9255
Shower primary prediction accuracy: 0.6000
Track fragment clustering accuracy: 0.9480
Interaction grouping accuracy: 0.9709
Particle ID accuracy: 0.8167
Primary particle score accuracy: 0.8864
e. Setup Evaluator¶
from analysis.classes.ui import FullChainEvaluator
# Only run this cell once!
evaluator = FullChainEvaluator(data, result, cfg, deghosting=True)
print(evaluator)
FullChainEvaluator(num_batches=10)
entry = 4 # Batch ID for current sample
print("Batch ID = ", evaluator.index[entry])
Batch ID = 4
3. Identifying Shower Primaries¶
Step 1: Get shower primary fragments¶
By using the primaries=True option, we can select out primary particles in this image. We will also load true_particles for comparison.
particles = evaluator.get_particles(entry, primaries=True)
true_particles = evaluator.get_true_particles(entry, primaries=True)
from pprint import pprint
pprint(particles)
[Particle( Batch=4 | ID=0 | Semantic_type: Shower Fragment | PID: Photon | Primary: 1 | Score = 98.09% | Interaction ID: 3 | Size: 1026 ),
Particle( Batch=4 | ID=2 | Semantic_type: Shower Fragment | PID: Photon | Primary: 1 | Score = 98.52% | Interaction ID: 3 | Size: 63 ),
Particle( Batch=4 | ID=3 | Semantic_type: Shower Fragment | PID: Electron | Primary: 1 | Score = 97.70% | Interaction ID: 3 | Size: 1509 )]
Alternatively, as you may have noticed, the primariness information is also stored in the Particle instance as an attribute with name is_primary. If you prefer to view the full image and then select out primaries manually:
particles = evaluator.get_particles(entry, primaries=False)
true_particles = evaluator.get_true_particles(entry, primaries=False)
from pprint import pprint
pprint(particles)
[Particle( Batch=4 | ID=0 | Semantic_type: Shower Fragment | PID: Photon | Primary: 1 | Score = 98.09% | Interaction ID: 3 | Size: 1026 ),
Particle( Batch=4 | ID=2 | Semantic_type: Shower Fragment | PID: Photon | Primary: 1 | Score = 98.52% | Interaction ID: 3 | Size: 63 ),
Particle( Batch=4 | ID=3 | Semantic_type: Shower Fragment | PID: Electron | Primary: 1 | Score = 97.70% | Interaction ID: 3 | Size: 1509 ),
Particle( Batch=4 | ID=4 | Semantic_type: Track | PID: Muon | Primary: 0 | Score = 67.55% | Interaction ID: 4 | Size: 1926 ),
Particle( Batch=4 | ID=5 | Semantic_type: Track | PID: Pion | Primary: 0 | Score = 34.94% | Interaction ID: 1 | Size: 197 ),
Particle( Batch=4 | ID=6 | Semantic_type: Delta Ray | PID: Photon | Primary: 0 | Score = 34.86% | Interaction ID: 4 | Size: 23 ),
Particle( Batch=4 | ID=7 | Semantic_type: Delta Ray | PID: Pion | Primary: 0 | Score = 84.00% | Interaction ID: 1 | Size: 25 ),
Particle( Batch=4 | ID=8 | Semantic_type: Delta Ray | PID: Proton | Primary: 0 | Score = 31.53% | Interaction ID: 4 | Size: 21 )]
Let’s quickly plot the particles and visualize which ones are predicted as primaries. Here is one way to do it with the trace_particles function:
from mlreco.visualization.plotly_layouts import white_layout, trace_particles, trace_interactions
traces = trace_particles(particles, color='is_primary', colorscale='rdylgn') # is_primary for coloring with respect to primary label
traces_true = trace_particles(true_particles, color='is_primary', colorscale='rdylgn')
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(traces, rows=[1] * len(traces), cols=[1] * len(traces))
fig.add_traces(traces_true, rows=[1] * len(traces_true), cols=[2] * len(traces_true))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
legend=dict(xanchor="left"),
autosize=True,
height=600,
width=1500,
margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)
The green voxels are predicted primary particles, while red indicates non-primary.
It is often easier to further break down the shower into different fragments and locate which of the shower fragments actually correspond to a predicted primary.
fragments = evaluator.get_fragments(entry)
Particle 12 has no PPN candidates!
Particle 13 has no PPN candidates!
Particle 18 has no PPN candidates!
Particle 22 has no PPN candidates!
Particle 23 has no PPN candidates!
Particle 29 has no PPN candidates!
Particle 31 has no PPN candidates!
Particle 37 has no PPN candidates!
Particle 40 has no PPN candidates!
Particle 43 has no PPN candidates!
Particle 45 has no PPN candidates!
traces = trace_particles(fragments, color='is_primary', colorscale='rdylgn') # is_primary for coloring with respect to primary label
traces_right = trace_particles(fragments, color='id', colorscale='rainbow') # This time, we'll plot the predicted particle
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(traces, rows=[1] * len(traces), cols=[1] * len(traces))
fig.add_traces(traces_right, rows=[1] * len(traces_right), cols=[2] * len(traces_right))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
legend=dict(xanchor="left"),
autosize=True,
height=600,
width=1500,
margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)
# TODO: Plot true fragment labels
Step 2: Identify the startpoint of the shower primary¶
During initialization of the Particle instance, PPN predictions are assigned to each particle if the distance between then is less than a predetermined threshold (attaching_threshold). PPN predictions that are matched to particles in this way are then stored in each Particle instance as attributes (ppn_candidates)
print("Minimum voxel distance required to assign ppn prediction to particle fragment = ", evaluator.attaching_threshold)
Minimum voxel distance required to assign ppn prediction to particle fragment = 2
fragments = evaluator.get_fragments(entry, primaries=False)
Particle 12 has no PPN candidates!
Particle 13 has no PPN candidates!
Particle 18 has no PPN candidates!
Particle 22 has no PPN candidates!
Particle 23 has no PPN candidates!
Particle 29 has no PPN candidates!
Particle 31 has no PPN candidates!
Particle 37 has no PPN candidates!
Particle 40 has no PPN candidates!
Particle 43 has no PPN candidates!
Particle 45 has no PPN candidates!
The first three columns are the \((x,y,z)\) coordinates of the PPN points. The fourth column is the PPN prediction score, and the last column indicates the predicted semantic type of the point.
We first visualize whether the predicted ppn candidates accurately locate the shower fragment start:
traces = trace_particles(fragments, color='id', size=1, scatter_ppn=True, highlight_primaries=True) # Set scatter_ppn=True for plotting PPN information
traces_true = trace_particles(true_particles, color='id', size=1)
fig = make_subplots(rows=1, cols=2,
specs=[[{'type': 'scatter3d'}, {'type': 'scatter3d'}]],
horizontal_spacing=0.05, vertical_spacing=0.04)
fig.add_traces(traces, rows=[1] * len(traces), cols=[1] * len(traces))
fig.add_traces(traces_true, rows=[1] * len(traces_true), cols=[2] * len(traces_true))
fig.layout = white_layout()
fig.update_layout(showlegend=False,
legend=dict(xanchor="left"),
autosize=True,
height=600,
width=1500,
margin=dict(r=20, l=20, b=20, t=20))
iplot(fig)
The left scatterplot highlighits primary shower fragments and its ppn candidates, while non-primaries are showed with faded color. The right plot shows true particle labels.
Identifying the primary shower fragments (as above) allow us to select all the voxels of the primary fragment which are close to the shower start, i.e. within some radius of the predicted PPN shower point. Of course, as expected from the scatterplot above, we may also include some cuts on the total voxel count to pick shower primary fragments that are large enough for our \(dQ/dx\) analysis.
For convenience, from now on we will only work with primary fragments:
fragments = evaluator.get_fragments(entry, primaries=True)
Particle 23 has no PPN candidates!
Step 3. Compute \(dQ/dx\) near the shower start¶
Let’s first fix some parameters for our \(dQ/dx\) computation. Let’s say the we select all points within a radius of 10 voxels from the predicted PPN shower start point of a given primary fragment and require that the selected segment size should at least be 3 voxels long.
from sklearn.decomposition import PCA
from scipy.spatial.distance import cdist
min_segment_size = 3 # in voxels
radius = 10 # in voxels
pca = PCA(n_components=2)
Write a compute_shower_dqdx function that takes a list of primary fragments and returns a list of computed dQ/dx values for each fragment.
def compute_shower_dqdx(frags, r=10, min_segment_size=3):
'''
Inputs:
- frags (list of ParticleFragments)
Returns:
- out: list of computed dQ/dx for each fragment
'''
out = []
for frag in frags:
assert frag.is_primary # Make sure restriction to primaries
if (frag.startpoint < 0).any():
continue
ppn_prediction = frag.startpoint
dist = cdist(frag.points, ppn_prediction.reshape(1, -1))
mask = dist.squeeze() < r
selected_points = frag.points[mask]
proj = pca.fit_transform(selected_points)
dx = proj[:, 0].max() - proj[:, 0].min()
if dx < min_segment_size:
continue
dq = np.sum(frag.depositions[mask])
out.append(dq / dx)
return out
compute_shower_dqdx(fragments)
[3761.2255977052623, 3888.1300989642573]
Step 4. Collect data over multiple images and plot results¶
iterations = 10
collect_dqdx = []
for iteration in range(iterations):
data, result = hs.trainer.forward(dataset)
evaluator = FullChainEvaluator(data, result, cfg, deghosting=True)
for entry, index in enumerate(evaluator.index):
# print("Batch ID: {}, Index: {}".format(entry, index))
fragments = evaluator.get_fragments(entry, primaries=True)
dqdx = compute_shower_dqdx(fragments, r=radius, min_segment_size=min_segment_size)
collect_dqdx.extend(dqdx)
collect_dqdx = np.array(collect_dqdx)
Segmentation Accuracy: 0.9851
PPN Accuracy: 0.8156
Clustering Accuracy: 0.9136
Shower fragment clustering accuracy: 0.9040
Shower primary prediction accuracy: 1.0000
Track fragment clustering accuracy: 0.9334
Interaction grouping accuracy: 0.9463
Particle ID accuracy: 0.8030
Primary particle score accuracy: 0.9153
Particle 13 has no PPN candidates!
Particle 42 has no PPN candidates!
Segmentation Accuracy: 0.9688
PPN Accuracy: 0.7979
Clustering Accuracy: 0.8609
Shower fragment clustering accuracy: 0.9102
Shower primary prediction accuracy: 0.7778
Track fragment clustering accuracy: 0.9569
Interaction grouping accuracy: 0.8940
Particle ID accuracy: 0.7564
Primary particle score accuracy: 0.8759
Particle 9 has no PPN candidates!
Particle 38 has no PPN candidates!
Particle 10 has no PPN candidates!
Particle 30 has no PPN candidates!
Particle 35 has no PPN candidates!
Segmentation Accuracy: 0.9932
PPN Accuracy: 0.8321
Clustering Accuracy: 0.9033
Shower fragment clustering accuracy: 0.9487
Shower primary prediction accuracy: 1.0000
Track fragment clustering accuracy: 0.9720
Interaction grouping accuracy: 0.9106
Particle ID accuracy: 0.8269
Primary particle score accuracy: 0.9273
Particle 29 has no PPN candidates!
Particle 35 has no PPN candidates!
Segmentation Accuracy: 0.9733
PPN Accuracy: 0.7615
Clustering Accuracy: 0.8377
Shower fragment clustering accuracy: 0.9434
Shower primary prediction accuracy: 0.9000
Track fragment clustering accuracy: 0.9534
Interaction grouping accuracy: 0.9414
Particle ID accuracy: 0.7778
Primary particle score accuracy: 0.9462
Particle 20 has no PPN candidates!
Particle 22 has no PPN candidates!
Particle 12 has no PPN candidates!
Particle 21 has no PPN candidates!
Segmentation Accuracy: 0.9814
PPN Accuracy: 0.7792
Clustering Accuracy: 0.9067
Shower fragment clustering accuracy: 0.9516
Shower primary prediction accuracy: 0.5455
Track fragment clustering accuracy: 0.9500
Interaction grouping accuracy: 0.9290
Particle ID accuracy: 0.8947
Primary particle score accuracy: 0.9154
Particle 12 has no PPN candidates!
Particle 14 has no PPN candidates!
Particle 49 has no PPN candidates!
Particle 51 has no PPN candidates!
Particle 16 has no PPN candidates!
Particle 30 has no PPN candidates!
Segmentation Accuracy: 0.9819
PPN Accuracy: 0.8271
Clustering Accuracy: 0.8859
Shower fragment clustering accuracy: 0.8787
Shower primary prediction accuracy: 1.0000
Track fragment clustering accuracy: 0.9530
Interaction grouping accuracy: 0.9193
Particle ID accuracy: 0.8868
Primary particle score accuracy: 0.9256
Particle 38 has no PPN candidates!
Particle 42 has no PPN candidates!
Particle 22 has no PPN candidates!
Particle 4 has no PPN candidates!
Particle 17 has no PPN candidates!
Particle 19 has no PPN candidates!
Segmentation Accuracy: 0.9840
PPN Accuracy: 0.8182
Clustering Accuracy: 0.8362
Shower fragment clustering accuracy: 0.8671
Shower primary prediction accuracy: 0.9333
Track fragment clustering accuracy: 0.9604
Interaction grouping accuracy: 0.9478
Particle ID accuracy: 0.9091
Primary particle score accuracy: 0.8868
Particle 4 has no PPN candidates!
Particle 19 has no PPN candidates!
Particle 19 has no PPN candidates!
Particle 10 has no PPN candidates!
Particle 14 has no PPN candidates!
Particle 36 has no PPN candidates!
Particle 39 has no PPN candidates!
Particle 17 has no PPN candidates!
Particle 25 has no PPN candidates!
Segmentation Accuracy: 0.9723
PPN Accuracy: 0.8146
Clustering Accuracy: 0.8628
Shower fragment clustering accuracy: 0.7941
Shower primary prediction accuracy: 1.0000
Track fragment clustering accuracy: 0.8976
Interaction grouping accuracy: 0.9419
Particle ID accuracy: 0.6957
Primary particle score accuracy: 0.8119
Particle 26 has no PPN candidates!
Particle 41 has no PPN candidates!
Particle 46 has no PPN candidates!
Segmentation Accuracy: 0.9494
PPN Accuracy: 0.8093
Clustering Accuracy: 0.8473
Shower fragment clustering accuracy: 0.8639
Shower primary prediction accuracy: 0.6364
Track fragment clustering accuracy: 0.9450
Interaction grouping accuracy: 0.9534
Particle ID accuracy: 0.8000
Primary particle score accuracy: 0.9008
Particle 10 has no PPN candidates!
Particle 16 has no PPN candidates!
Particle 22 has no PPN candidates!
Segmentation Accuracy: 1.0000
PPN Accuracy: 0.6214
Clustering Accuracy: 0.9372
Shower fragment clustering accuracy: 1.0000
Shower primary prediction accuracy: 0.0000
Track fragment clustering accuracy: 1.0000
Interaction grouping accuracy: 0.9000
Particle ID accuracy: 1.0000
Primary particle score accuracy: 0.8000
collect_dqdx
array([ 8383.44787459, 3075.15671474, 10275.39597007, 5759.81818464,
2048.33477851, 6483.26586487, 6711.89453074, 4442.57608566,
5810.30288252, 2084.21310505, 15320.11624131, 5330.756099 ,
4077.75588921, 3100.59967366, 3975.09149345, 3203.11402457,
4501.66203046, 2539.68349723, 4118.45226987, 5259.90855816,
4267.68260069, 3879.05482131, 5105.65571194, 5702.04481627,
3774.81705912, 7069.72264419, 3782.59048981, 2857.099469 ,
4209.95435594, 2411.26060632, 1770.08554637, 3933.01448637,
3560.54564803, 5301.89095261, 4443.64398492, 5435.98994269,
3248.57105291, 3641.49973808, 4909.68242799, 1937.21630116,
3668.666896 , 2377.0897961 , 5423.19395053, 3612.28862708,
1348.58141929, 2295.38587576, 2251.51456024, 3672.65459286,
3254.90768849, 4100.30846616, 4982.4629635 , 5855.32706474,
4145.05137581, 4346.63280712, 11732.50798136, 6856.72447656,
3852.3753254 , 2342.01964364, 4185.40556866, 4929.1805666 ,
1585.08337361, 6748.9593629 , 1394.22589888, 1473.70338119,
7568.51794417, 4242.0805526 , 1660.48735056, 1260.2668119 ,
3858.60545567, 3665.90093946, 4916.78392715, 4581.43252763,
1997.71039713, 4878.08911417, 5269.80159383, 7909.49298996,
2716.08654783, 7782.43046032, 4730.61239327, 2540.66956857,
2158.67235926, 4567.52157391, 1908.14601307, 4567.61153747,
8685.72533064, 5874.09140703, 5191.60880059, 3156.7778149 ,
1644.08159374, 2898.62098508, 4306.86633509, 4433.66329358,
2065.59814292, 4443.62302853, 2768.49099097, 7208.52343767,
3542.01472489, 1349.81173533, 2347.78723609, 1615.27970814,
4397.25836733, 4680.11439817, 6320.96339051, 5705.68633357,
2839.11812281, 2724.78354291, 2481.55748291, 4066.90042351,
2178.35369663, 2626.39918086, 1706.91480243, 3277.94537695,
6610.11741142, 2897.46137901, 2124.82437166, 1908.56779833,
2397.91343631, 3105.56355578, 3793.02011089, 2374.60806139,
5266.3662759 , 7004.7296959 , 5824.07798361, 2037.67489369,
7214.89858497, 1794.39217209, 11365.73869412, 3990.37565522,
5523.23893662, 5077.7209911 , 2830.31733061, 6690.00966216,
2277.16182534, 2123.06629568, 1906.32678771, 2045.40149631,
4825.26648769, 3426.94249517, 3139.10820024, 1722.74570783,
4668.92211927])
import matplotlib.pyplot as plt
import seaborn
seaborn.set(rc={
'figure.figsize':(15, 10),
})
seaborn.set_context('talk')
plt.hist(collect_dqdx, range=[0, 10000], bins=50)
plt.xlabel("dQ/dx")
plt.ylabel("Predicted primary shower fragments")
Text(0, 0.5, 'Predicted primary shower fragments')
